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1 Introduction 



In the 1980s, Arak and Surgailis introduced a class of planar Markov fields whose 
realisations form a coloured tessellation of the plane. The basic idea is to use the 
lines of an isotropic Poisson line process as a skeleton on which to draw polygonal 
contours with the restriction that each line cannot be used more than once. Note 
that many tessellations can be drawn on the same skeleton and that the contours may 
be nested. The polygons are then coloured randomly subject to the constraint that 
adjacent ones must have different colours. Formally, the probability distribution of such 
polygonal Markov fields is defined in terms of a Hamiltonian with respect to the law 
of the underlying Poisson line process, which can be chosen in such a way that many 
of the basic properties of the Poisson line process (including consistency, Poisson line 
transects, and an explicit expression for its probability distribution on the hitting set 
of bounded domains) carry over and a Markov property holds. Another useful feature 
is that a dynamic representation in terms of a particle system is available. See [21 13] 
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for further details and [1] or [27] for alternative point rather than line based models. 
The special case where all vertices have degree three was studied in pE]. 

The simplest and most widely studied example of a planar Markov field is the Arak 
process pj which consists of self-avoiding closed polygonal contours. Hence interaction 
is restricted to a hard core condition on the contours, and there are exactly two colour- 
ings using the labels 'black' and 'white' such that adjacent polygons have different 
colours. For the Arak process, the Hamiltonian is proportional to the total contour 
length by a factor two. One may introduce further length-interaction by changing the 
proportionality constant. Doing so, Nicholls [IB] and Schreiber [2T] consider the separa- 
tion between the black and white regions as the interaction gets stronger; Van Lieshout 
and Schreiber [H] develop perfect simulation algorithms for these models. For the more 
general case in which both length and area terms feature in the Hamiltonian, Schreiber 
[20] proposes a Metropolis-Hastings scheme based on the dynamic representation of 
the Arak process, which Kluszczyhski et al. [11] adapt and implement to solve fore- 
ground/background image segmentation problems. An anisotropic Arak process can 
be defined through local activity functions instead of the length functional [221 121] and 
allows for increased flexibility while preserving desirable basic properties including a 
dynamic representation. This representation forms the basis of a stochastic optimisa- 
tion algorithm in the context of image segmentation, implemented by Matuszak and 
Schreiber [17], and helps to gain insight into the higher order correlation structure [22] . 

Polygonal Markov field models with contours that may also be joined by a vertex of 
degree three or four [21 [3] are much less well-understood due to the more complicated 
interaction structure. Papers in this direction include [71 [T21 [T9] . 

In a previous paper [23], we introduced a class of binary random fields that can be 
understood as discrete counterparts of the two-colour Arak process. The aim of the 
present paper is to extend this construction to allow for an arbitrary number of colours 
and to relax the assumption in [23] that no polygons of the same colour can be joined by 
corners only. Our construction is two-staged: first a collection of lines is fixed to serve as 
a skeleton for drawing polygonal contours (a regular lattice being the generic example), 
then the resulting polygons are coloured in such a way that adjacent ones do not have 
the same colour. The analogy with continuum polygonal Markov fields is exploited 
to define Hamiltonians that are such that the desirable properties of these processes 
mentioned above hold. Moreover, we propose new simulations techniques that combine 
global changes with the usual local update methods employed for random fields on finite 
graphs |29j. It should be stressed, though, that the discrete models considered in this 
paper are not versions of continuum polygonal Markov fields conditioned on having 
their edges fall along a given finite collection of lines. 

The plan of this paper is as follows. In Section [2] we define a family of admissible 
multi-colour polygonal configurations built on regular linear tessellations, and define 
discrete polygonal fields with special attention to consistent ones. In Section [3] we 
present a dynamic representation of consistent polygonal fields, which is used to prove 
the main properties of such models. In Section HI we exploit the dynamic represen- 
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tation to develop a simulation method for consistent polygonal fields. The method 
is generalised to arbitrary Gibbs fields with polygonal realisations and applied to the 
detection of linear networks in images in Section |5l We conclude with a discussion. 

2 Random fields with polygonal realisations 

First, recall the definition of a regular linear tessellation from [23] . 

Definition 1. A regular linear tessellation of the plane is a countable family T of 
straight lines in such that no three lines of T intersect at one point and such that 
any bounded subset of the plane is hit by at most a finite number of lines from T ■ 

For a bounded open convex set T induces a partition of D into a finite collection 
Dj- of convex cells of polygonal shapes, possibly chopped off by the boundary. Below 
we shall always assume that the boundary dD of D contains no intersection points of 
lines from T and that the intersection of each line / G T with dD consists of exactly 
two points. To each line /, we ascribe a fixed activity parameter vr^ G (0, 1) to allow for 
the possibility to favour some lines over others. 

The next step is to assign a colour to each of the convex cells in the partition of D 
induced by the lines in T. Write {1, . . . , fc} for the set of colour labels. In this paper, we 
concentrate on the case that k > 2. Such a colouring gives rise to a graph whose edges 
are formed by the boundaries between cells that have been assigned different colours. 
For technical convenience, we shall assume that edges are open, that is, they do not 
contain the vertices in which they intersect. Faces of the graph, which are unions of 
cells of Df, are said to be adjacent if they share a common edge. 

Definition 2. The family f^(T) o/ admissible coloured polygonal configurations in D 
built on T consists of all coloured planar graphs 7 in the topological closure D = DUdD 
of D such that 

• all edges lie on the lines ofT; 

• all interior vertices, i.e. those lying in D, are of degree 2, 3 or 4; 

• all boundary vertices, i.e. those lying on dD, are of degree 1; 

• no adjacent faces share the same colour. 

Throughout this paper, the notation 7 is used for (admissible) planar graphs and the 
hat notation 7 for the graph with colours assigned to its faces. In this notation, VoiT) 
stands for the family of all planar graphs 7 arising as interfaces between differently 
coloured faces in 7 G r/3(T). Note that for the case k = 2 treated in [23], all interior 
vertices have degree two. Vertices of degree two are also known as V-vertices, those of 
degree three as T-vertices, and vertices of degree four as X-vertices. 

To avoid confusion, it is important to distinguish between D-j- and the members 
of r£)(T), even though they all partition D. In the sequel we shall reserve the terms 



3 



> > 



Figure 1: Admissible polygonal configuration. 



'segments' and 'nodes' for the former, 'edges' and 'vertices' for the latter. Thus, nodes 
are intersection points of lines in T, which are joined by segments. Edges of 7 are 
maximal unions of connected coUinear segments that are not broken by other such 
segments lying on the graph (corresponding to vertices of degree three or four). Pri- 
mary edges are maximal unions of connected coUinear segments, possibly consisting of 
multiple edges due to T- or X-vertices. Likewise, a vertex of 7 is a point where two 
edges of 7 meet or where an edge of 7 meets the boundary dD, nodes lying on the 
interior of graph edges are not considered to be vertices. These concepts are illustrated 
in Figure [H The family T consists of two orthogonal line bundles and induces diamond 
shaped cells. One element oiVniT) is plotted in Figured! Consider the polygonal face 
indicated by 'C that is chopped off by the boundary. There are 2 boundary vertices 
and 12 interior vertices, 5 of degree 2, 7 of degree 3 and none of degree 4. We count 
17 nodes on 16 complete segments, with 2 segments partly visible due to truncation 
by the boundary. C has 13 edges and its boundary lies on 8 primary edges. Note that 
in the literature, primary edges are sometimes also known as sides. See e.g. pS] for a 
discussion on nomenclature for tessellations. 

We are now ready to define (discrete) polygonal fields. 

Definition 3. Let T he a regular linear tessellation and ascribe fixed activity param- 
eters Til G (0, 1) to the lines I G T. For a function Hd : ^^(T) t-)- M U {+00}, the 
(discrete) polygonal field Ay,o with Hamiltonian T-Ld is the random element in roiT) 
such that 

Here E*{'y) denotes the set of primary edges in 7 and l[e] G T is the straight line 
containing the open edge e. 

The constant 

mD]= Yl ^M-'Hom n ^iie] (2) 

eetoiT) e&E*{e) 



that ensures that P is a probabihty distribution is called the partition function. Note 
that the polygonal field is a Gibbs field with Hamiltonian 

noil)- 

ee-E*(7) 

The terms log 7ii[e] represent the energy needed to create the edges. We prefer to use the 
term polygonal field since the consistent fields to be considered shortly are inspired by 
the Arak-Surgailis polygonal Markov fields in the continuum, and, more importantly, 
because the graph-theoretical formulation leads us to define novel simulation techniques 
for discrete random fields. 

A careful choice of Hamiltonian in ([T]) leads to consistent polygonal fields. Recall 
that k is the number of colour labels. 

Definition 4. Let ay G [0, 1] and set 

I k-2 \ ay k-2 

The polygonal fields (J\) defined by Hamiltonians of the form 

$d(7) = -Ny{-i)\ogay-NT{i)\og{{k-l)aT)-Nx{l)\og{{k-l)ax) 
+ card(E(7)) log(A; - 1) 

-EE Ml-evr,)+ E (l - ^vr.^Tr,,] (3) 

are referred to as consistent polygonal fields. Here, Ny, Nt and Nx denote the number 
of V-, T- and X-vertices, caid{E{j)) is the number of edges in 7, n(/i,/2) G 7 ranges 
through the nodes of T that either lie on edges of 7 or coincide with one of its vertices, 
and I ^ e means that the line I intersects but is not collinear with e. We use the 
convention that x 00 = 0. 

A few remarks are in order. The parameter ay controls the relative frequency of 
V^-vertices, ar and ax that of vertices of degrees three and four respectively. These 
parameters are not independent and, given the number of colours k, ax and ax are 
uniquely determined by ay. This dependence will become more explicit in the dynamic 
representation to be derived in Section [31 Typical realisations for = 3 and ay equal 
to 0, 1/2 and 1 are shown in Figure |2i In the left-most panel, ay = and there no 
vertices having degree two; in the right-most panel, ax = so that no vertices have 
degree four. The central panel displays a coloured configuration with vertices of all 
degrees. Visually, the three patterns are strikingly different. 

For now, let us consider the special case k = 2 and ay = 1. First note that the 
family r£)(T) of admissible coloured polygonal configurations does not include any 
member with an interior vertex of degree three. Therefore = almost surely. 
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Figure 2: Realisations of »4$^ with ay = (left), ay = 1/2 (centre) and ay = 1 (right). In 
all panels, the number of colour labels k = 3. 

Moreover, as ax = 1 — ay = 0, the probability of any 7 that contains X-vertices is 
zero. Hence, almost surely all interior vertices have degree two and E*{'j) = -£^(7). 
Moreover, \og{k — 1) = log ay = 0, and simplifies to 

cf. [23]. 

The nomenclature in Definition |4] is justified by the following result. 

Theorem 1. The polygonal field A^^^ is consistent: For hounded open convex D' C 
C M^, the field A^j^ fl D' coincides in distribution with Ai^,^,. 

By letting D increase to M^, Kolmogorov's theorem allows us to construct the whole 
plane extension of the process A^ such that the distribution of coincides with that 
of Aq, n D for all bounded open convex D <ZM?. The proof of Theorem [1] relies on a 
dynamic representation, which is the topic of the next section. 

3 Dynamic representation of consistent polygonal fields 

Below, we present a dynamic representation for discrete consistent polygonal fields 
in analogy with the corresponding representation in [21 Sections 4 and 5]. The idea 
underlying this construction is to represent the edges of the polygonal field as the 
trajectory of a one-dimensional particle system evolving in time. More specifically, we 
interpret D as a set of time-space points (t, y) E D and refer to t as the time coordinate, 
to y as the spatial coordinate of a particle. In this language, a straight line segment in 
D stands for a piece of the time-space trajectory of a moving particle. Compared to 
the bi-coloured case discussed in [23], note that in the current multi-colour context, in 
general it is not sufficient to consider colourless graphs only and assign colourings with 
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equal probability afterwards. Moreover, the interaction structure is not restricted to a 
hard core condition. 

For convenience, we assume that no line in T or segment of dD is parallel to the 
spatial axis. Since we might simply rotate the coordinate system otherwise, these 
assumptions do not lead to a loss of generality. Recall that by assumption, each line 
/ of T intersects the boundary at two points. The two intersection points are ordered 
according to time and the one with the smaller time coordinate is denoted in(/,Z)). 
Furthermore, no three lines of a regular linear tessellation intersect in a single point. 

To define the dynamics, the left-most point of D is assigned a random colour chosen 
uniformly from the k possibilities. Let particles be born independently of other particles 

• with probability ayTfi-i^iTi^/lk — 1) at each node n{li,l2), that is, the intersection 
of two lines li and I2 in T, which falls in D (interior birth site), 

• with probability 7r;/(l + tt/) at each entry point in(/,D) of a line / G T into D 
(boundary birth site). 

Each interior birth site n{li,l2) emits two particles moving with initial velocities such 
that the initial segments of their trajectories lie on the lines li and I2 of the tessellation 
that emanate from the birth site, unless another particle (either a single one or two 
colliding particles) previously born hits the site, in which case the birth does not occur. 
Each boundary birth site in(/, D) emits a single particle moving with initial velocity 
such that the initial segment of its trajectory lies on /. Note that no precaution similar 
to the one for interior birth sites above is needed because boundary birth sites cannot 
be hit by previously born particles. The initial trajectory or trajectories of a birth 
event bound a new polygonal region, the colour of which is chosen uniformly from the 
k — 1 colours that differ from that of the polygon just prior to the birth, or in other 
words, lying to the left of the birth site. 

All particles evolve independently in time according to the following rules. 

(El) Between the critical moments listed below each particle moves freely with con- 
stant velocity. 

(E2) When a particle hits the boundary dD, it dies. 

(E3) In case of a collision of two particles, that is, equal spatial coordinates y at some 
time t with {t, y) = n{li, Ij) G D, distinguish the following cases. 

a If the colours above and below (t, y) are identical, say i, with probability ay 
both particles die. With probability ax both particles survive to create a 
new polygon whose colour is chosen uniformly from those not equal to i (cf. 
Figure [3]) . 

b If the colours above and below (t, y) are different, say i and j, with probability 
ot, each of the two particles survives while the other dies. With probability 
[k — 2)axl{k — 1) = 1 — 2aT, both particles survive to create a new polygon 
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Figure 3: Top row ((E3) [a]): In the left-most panel, both particles die; they survive in the 
panel on the right and create a new polygon coloured blue. Bottom row ((eE3) [b]): In the 
left-most panel, one particles survives; both particles survive in the panel on the right and 
create a new polygon coloured red. 

whose colour is chosen uniformly from those not equal to either i or j (cf. 
Figure [3]) . 

Recall that a collision prevents a birth from happening at the node. 

(E4) Whenever a particle moving in time-space along k E T reaches a node n{li, Ij), it 
changes its velocity so as to move along Ij with probability ayT^iJik — 1), it splits 
into two particles moving along /j and Ij with probability {k — 2)aTT^iJ{k — 1) 
and keeps moving along /j otherwise (with probability 1 — evr;^.). In case of a 
split, a new polygon is created whose colour is chosen uniformly from the k — 2 
possibilities. See Figure HI 

The dynamics described above define a random coloured polygonal configuration 
T>£). The key observation is that its distribution is identical to that of Aq,^. 

Theorem 2. The random elements and coincide in distribution. 

Proof: In order to calculate the probability that some 7 G T is generated by 
the particle dynamics El— E4, observe that 

• each edge e G -£^(7) whose initial (lower time coordinate) vertex lies on dD yields 
a factor 7r;[e]/(l + 7r;[e]) (boundary birth site) times Hier i^oe^^ " ('^^ veloc- 
ity/split updates along e) times l/(fc — 1) for the colour; 

• the two edges ei, 62 G £'(7) emanating from a common interior birth site n(/i, I2) 
yield a factor avT^hT^i2/ {k — 1) (the birth probability) times IlLi Y\ieT (-'-"^^0 
(no velocity/split updates along Cj) times l/{k — 1) for the colour; 
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Figure 4: (E4): Velocity update (left), split (middle) and continuation (right) of a trajectory 
outlined in red. 

• each edge e G -£^(7) arising due to a velocity update yields a factor Q;y7r;[e]/(/c — 1) 
(velocity update probability) times YlieT i^e(-'- ~ (^^ velocity/split updates 
along e); 

• the two edges 61,62 G -£'(7) arising from a split at node nili^l^) where 61 is a 
continuation and 62 a new direction yield a factor {k — 2)aT7^i2/{k — 1) (the split 
probability) times IlLi IliGr I'^eS^ " ^'^'-^ ^^'^ velocity/split updates along 6^) 
times l/{k — 2) for the colour; 

• the two edges 61,62 G -£^(7) arising due to an X-coUision contribute a factor 
dx/ik — 1) times IlLi Y\i^t I'^^eS^ ^ (^^ velocity/split updates along 6j); 

• each edge e emanating from a T-coUision contributes a factor ar times Ylier i'>oei^~ 
eni) (no velocity/split updates along e); 

• each collision of V-type contributes a factor ay; 

• the initial colour choice contributes a factor l/k; 

• the absence of birth sites in nodes n(/i, ^2) £ that do not belong to 7 yields the 
factor nn(Zi,Z2)GZ3\7(l - ayvr;,7ri2/(/i; - 1)) ; 

• the absence of boundary birth sites at those entry points into D of lines of T which 
do not give rise to an edge of 7 yields the factor Hier, «nD^0, m(«,D)07(l + 

Collecting all factors implies that the probability of 7 is 

(n no--))( n 0-^-.-.))x( n iX^) 

times 

which completes the proof. □ 
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As an immediate consequence of the proof of Theorem [2l we obtain an exphcit and 
simple expression for the partition function of consistent polygonal fields. 

Corollary 1. The partition function ^ is given by 

n n (i-fc^^'^^O ' 

Having established a proper dynamic representation, Theorem [1] follows from The- 
orem [2] in complete analogy with the proof of [21 Thm. 5.1] as in [23, Thm. 1]. A 
discrete analogue of the Poisson line transect property Thm. S.l.c] holds as well. 
Indeed, combining consistency with the boundary birth mechanism of the dynamical 
representation, we obtain the following. 

Corollary 2. Letl he a straight line that contains no nodes ofT ■ Then, the intersection 
points and intersection directions of I with the edges of the polygonal field A<^ coincide 
in distribution with the intersection points and directions of I with the line field Aj- 
defined to be the random sub- collection of T where each straight line I* E T is chosen 
to belong to Aj- with probability tti* / (l+vrp) and rejected otherwise, and all these choices 
are made independently. 

To conclude this section, we turn to Markov properties that are direct consequences 
of the dynamic representation. 

A spatial Markov property reminiscent of that enjoyed by the Arak-Clifford-Surgailis 
model in the continuum is the following. For a piece wise smooth simple closed curve 
^ C containing no nodes of T, the conditional distribution of A<^ in the interior of 
9 depends on the configuration exterior to 9 only through the intersections of 9 with 
the edges of the polygonal field and through the colouring of the field along 9. 

To relate our model to Gibbs fields commonly used in image analysis, assume for 
the remainder of this section that T forms a regular lattice and D is an m x n rectangle. 
In this case, D is divided by T in square cells, known as pixels. There is a one-to- 
one correspondence between a coloured polygonal configuration 7 and the array of 
pixel colours. Indeed, the colour of a pixel is that of the face of 7 that it falls in, 
and, reversely, the edges of 7 are composed of the segments between pixels of different 
colours. In this dual framework, we obtain the following local Markov factorisation. 

Corollary 3. Let D be an m ^ n array and let T = {/i, . . . , /m+n-2} be the cor- 
responding regular linear tessellation with the indices chosen in such a way that for 
i = 1, ... ,n — 1, li is the horizontal line between the i-th and {i + l)-st row and for 
i = 1, . . . ,m — 1, In+i-i is the vertical line between the i-th and {i + l)-st column. 
Then the random element /1$^ is the dual of a random vector X = {Xi, . . . , Xmn) 
of pixel values indexed in column major order with a joint probability mass function 
P(7) = F{Xi = xi] ■ ■ ■ ; Xmn = Xmn) that factorises as 

n m—1 

= xi) ]^P(Xj = Xi I Xj_i = JJ^ P(Xj„_|_i = Xin+i I X{i_i)n+i = a;(i_i)„+i) 

i=2 i=l 
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m— 1 n 

for Xi E {1, . . . , k}, i = 1, . . . , mn. 

Proof: Choose the time direction in the dynamic representation in such a way that 
the chronological order of the nodes coincides with the column major order and recall 
that the conditional behaviour at a node depends only on the colours and trajectories 
immediately to its 'left', i.e. immediately prior to it in time. For the first pixel, = 
Xi) is the probability that the initial colour is Xi, which is 1/A; as this colour is chosen 
uniformly from the set {1, . . . ,k}. Next, the probabilities of the pixel values in the first 
column are derived from the boundary birth mechanism. Thus, for Xi, Xi^i e {1, . . . , A;} 
and i — 2, . . . ,n, 



1 '^',-1 



F{Xi = Xi I = Xi-i) = 



if Xi ^ 



k~l l+7ri._j 



Similarly, for Xin+i, e {l,...,k} and i = 1, . . . , m - 1, 

I l+TT, — ^(i— l)n+l- 

V 'n+i— 1 

Use the shorthand notation Fij{x \ u,v,w) for 

Then, ior u ^ v w u, 

TO ^ I \ \ <^v ilx^v 

^■'^ ' ' ' ^ I axlik-^) \ix^v 

by (E3a), and 

r,,Kx\u,v,w)-^^^l^^_^s^ ifx^{T;,w} 
by (E3b) regardless of i^j. Furthermore, the interior birth mechanism implies that 

TO / I \ — \ ctv"^/ -i^i + ^ 1)^ \i X ^ u 

Vi^yx I u,u,u)-^^_ ^^^^l^^^^^^j - 1) iix = u 

and finally (E4) determines the probabilities 

1 — eTx^_^ ii X — V 

Pjj(a; \u,u,v) — avTi^_^/{k — 1) if a; = -u, 

ar7r/^_i/(A) - 1) \i x ^ {u,v} 
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,_i ifx = v; 

/{k — 1) if X = u] 
/{k-l) if X ^ {u,v]. 



Collecting all terms completes the proof. 



□ 



Models for which the above factorisation holds have been dubbed mutually compat- 
ible Gibbs random fields by Goutsias [9J. In particular, the interior birth mechanism, as 
well as the collisions and path propagation described in the dynamical representation's 
(E3)— (E4) correspond to the local transfer function in [9], see also[5]. 

Note that we have chosen the time direction so as to conform to column major 
order. A fortiori, such models are Markov random fields with the second order neigh- 
bourhood structure in which horizontally, vertically, and diagonally adjacent pixels are 
neighbours. A similar factorisation holds for any choice of the time axis. Moreover, 
there is no need for D to be a rectangle. Indeed, because of the assumptions on T, every 
interior node is hit by exactly two segments that are adjacent to three, not necessarily 
rectangular, pixels. See Figure [H] for an example. The notation, however, becomes 
more cumbersome. For this reason, we prefer to use the graph-theoretical formulation 
with its neater formulae. 

4 Birth and death dynamics 

In this section, we use the dynamic representation of Section [3]to propose dynamics that 
are reversible and leave the law of invariant. These dynamics will serve as stepping 
stone for building Metropolis-Hastings dynamics for the general polygonal field models 
([T]). For the two-colour case, algorithms inspired by dynamic representations can be 
found in [HI 1201 ESI IIZ]- In that case, however, it is sufficient to focus on colour blind 
polygonal configurations as the colouring is completely determined by the colour at a 
single point. For k > 2, this is no longer the case and we have to explicitly incorporate 
colours in our dynamics. Moreover, particles do not necessarily die upon collision and 
the disagreement loop principle of Schreiber [20] no longer applies. 

The basic continuous time dynamics we propose consist of adding and deleting 
particle birth sites. In order to fully explore the state space, recolouring will also be 
necessary, at some fixed rate r > 0. We work with a constant death rate 1. The birth 
rate at a boundary entry point in{l,D) and a vacant internal node n{li,l2) are set to 
7ii and ay7r/^7r;2/(fc — 1 — ayTT/^vr^J to satisfy the detailed balance equations 



respectively 



k-l 



l + Tll 



T^hT^h X 1 



X 1 




TT/jTTjj J X birth rate(n(/i, /2))- 



) 



X birth rate(in(Z, D)) 
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Recall that if nili^h) is hit by some previously bom particle, the birth is discarded. 
For computational convenience, we shall keep track of the discarded births during the 
dynamics. 

In case of a birth update, the particle (s) emitted by the birth site are given tra- 
jectories in accordance with E4. Upon collisions, E3 is invoked. Whenever possible, 
existing trajectories are re-used. A dual reasoning is applied to deaths. 

To make the above ideas precise, suppose that the current state is 7, understood 
here to include the knowledge of all discarded birth sites, which we modify by adding 
or deleting a (discarded) birth site to obtain 7'. We shall use the following segment 
classification: 

plus the segment does not lie on any edge of 7 but it does lie on some edge of 7'; 

minus the segment lies on some edge of 7 but it does not lie on any edge of 7'; 

changed the segment lies on some common edge of 7 and 7' but the colour of at least 
one of its adjacent polygonal faces has changed, or the segment lies in the interior 
of faces of 7 and 7' having different colours. 

The dynamics are now as follows. In case of a birth at node n{li,l2) with coor- 
dinates {t,y), two plus segments arise along li and I2 forward in time. In case of a 
boundary birth, a single plus segment is generated. Similarly, in case of a death, one 
or two minus segments arise. We order the nodes with first coordinate larger than t in 
chronological order and update them one at a time until some further time t' > t for 
which the intersections of 7 and 7' with the vertical line specified by first coordinate t' 
are identical. 

At each node, we first check whether the node is hit either by some segment marked 
'plus' or 'minus', or by some 'changed' segment of 7. When using the word 'hit' we 
shall always mean that the tail of the segment has a smaller time coordinate than the 
node at its head. We shall use the phrase 'emanating segments' for those segments 
whose tail is at the node and whose head has a larger time coordinate than the node. 
If no such segment exists, we need to check whether the node is a non-discarded birth 
site. If it is, due to e.g. nesting, the colour just prior to the birth site may be different 
in 7 and 7'. In this case, a new colour is chosen for the region to the right of the node 
from those not equal to the colour just prior to the node. Otherwise, the status quo is 
propagated. 

If the node is hit by two marked segments ('plus', 'minus' or 'changed' and belonging 
to 7), or by one such segment and one that is an unchanged common edge of 7 and 7' 
(which we label 'old') a collision update is made as outlined in Table [TJ 

In the remaining case that the node is hit by a single segment marked either 'plus' 
or 'minus' or by a segment of 7 labelled 'changed', the path is updated as outlined in 
Table H 

An illustration is given in Figure O The current node n to be updated is in the 
middle of the panels. In 7, the node is a birth site: no segments hit n and there are 
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minus/minus 
at vacant node 


label emanating segments; 


minus / minus 

at discarded birth site 


implement birth by choosing new colour from those 
not equal to that prior to the discarded birth site; 
label emanating segments; 


minus/old 
minus/plus 
minus / changed 


invoke E4; 

label emanating segments; 


plus/old 
plus/changed 
old/changed 
plus/plus at 
vacant node 


invoke E3; 

label emanating segments; 


plus/plus 
at birth site 


discard the birth; 

invoke E3 and label emanating segments; 


changed / changed 


check whether the colours above and below the face 
bounded by the two hitting segments agree in 7 and 7'; 
if so, do nothing; 

otherwise invoke E3 Mud \i\hc\ (nnaualiug segments. 


Table 1: CoUision updates. 



plus path 
at birth site 


discard the birth; 
invoke E4; 

label emanating segments; 


plus path 
at vacant site 


invoke E4; 

label emanating segments; 


minus path 

at discarded birth site 


implement birth by choosing new colour from those 
not equal to that prior to the discarded birth site; 

label emanating segments; 


minus path 
at vacant site 


label emanating segments; 


changed path 


in case the path splits, choose a new colour 

from those not equal to the colours above and below the path; 

label emanating segments. 



Table 2: Path updates. 
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two emanating segments. In the new polygonal configuration 7', the node is hit by 
a segment separating the green from the blue face which is therefore labelled 'plus'. 
According to Table [21 in 7', the birth gets discarded and we invoke E4, say resulting 
in the decision to split. Hence, we must also choose a colour other than green or blue 
for the region to the right, that is, forwards in time, for example red. The labels of the 
two emanating segments are 'changed' for the one separating the blue and red faces, 
and 'old' for the one forming the boundary between the red and green faces. 




Figure 5: Left: 7. Right: 7'. 

For recolouring, classic local colour switches are used, as detailed in for example 
[29] . As for the two-colour case in [23], we obtain the following result. 

Theorem 3. The distribution of the consistent polygonal field A<i>jy is the unique invari- 
ant probability distribution of the birth- death-recolour dynamics described above upon 
ignoring discarded birth attempts, to which they converge in total variation from any 
initial state 7 G rz)(T) for which F{A^^ =7) > 0. 

Proof: The case k = 2 was considered in [23]. Hence assume k > 2. The total 
transition rate is bounded from above by a positive constant. Indeed, for each internal 
birth site n(/i,/2), the rate 

avTihTTiJ{k - 1) ^ l/{k - 1) ^ 1 ^ ^ 
I- ayKi^TTiJik-l) - 1-1/ {k-l) A; - 2 " 

The death rate for a (discarded) birth at n(Zi, I2) is equal to 1. Similarly, the birth rate 
TTi at entry points in(l, D) = tti < 1 and the death rate is again equal to 1. Therefore, 
an upper bound is the sum of r, the number of nodes and the number of lines hitting 
D. Thus, our dynamics can be algorithmically generated by a Poisson clock of con- 
stant rate, and an embedded Markov transition matrix that governs the transitions. 
This transition matrix restricted to admissible coloured polygonal configurations with 
positive probability under the putative invariant probability distribution is irreducible, 
since any state can be reached from any other state by successively removing all birth 
attempts, choosing an appropriate initial colour and then building the target state by 
successively adding particles. Hence the dynamics constitute a finite state space ir- 
reducible Markov process and there exists a unique invariant probability distribution. 
See for example Theorem 20.1 in [13]. The same theorem also yields the converge in 
total variation. The invariance of the distribution including discarded birth attempts 
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follows from the invariance of the Bernoulli birth site probabilities under the dynamics, 
the fact that the trajectories preserve El— E4 by design, and the well-known invari- 
ance of the local colour switches [29( Section 10.2]. Summing over any discarded birth 
attempts completes the proof. □ 

In fact, the dynamics are reversible. Consider for example the collision of a positive 
path with a birth site as in Figure El In this example, the transition rate is multiplied by 
{k — 2)aT'n'i./{k — l) for the split and l/{k — 2) for the colour, amounting to ax'n'i^/ (k — l) 
whereas the probability of 7 including the knowledge of discarded births gains a factor 
[k — l)aT'n'ij /{k — 1) = utt^Ij (the first term A; — 1 to compensate for the choice of colour 
in 7, the second for the split and colour in 7'). The reverse collision of a minus path 
with a discarded birth site yields a factor l/(fc — 1) for the rate. Thus, this type of 
collision is reversible. Similar calculations can be made for the other types of collision 
and are left to the reader. 

For general polygonal field models ([T]) with a Hamiltonian that is the sum of ([3]) 
and some other term T/^i, we propose a Metropolis-Hastings algorithm. The algorithm 
has the same birth, death and recolour rates as the dynamics presented in Section HI 
The difference is that a new state 7' is accepted with probability 

min(l,exp[HB(7)-Hzp(7')]), (4) 

whereas the old state is kept with the complementary probability. An example oiT-Lu 
will be presented in Section O 

Theorem 4. Let T-Ld '■ ^d{T) h- M &e finite. Then the distribution of the polygonal 
field A^jj+Ud unique invariant probability distribution of the Metropolis-Hastings 

dynamics. Moreover, these dynamics converge in total variation to from any 

initial state 7 G Td(T) for which P(^$^ =7) > 0. 

Proof: By the assumption on the Hamiltonian and arguments analogous to those 
in the proof of Theorem [3l the embedded Markov transition matrix that governs the 
transitions is irreducible. Hence the dynamics constitute a finite state space irre- 
ducible Markov process and there exists a unique invariant probability distribution, cf. 
Theorem 20.1 in [13], to which they converge in total variation. By Theorem |31 the 
birth-death- recolour dynamics leave the distribution of Aq,j^ invariant. The modifica- 
tion by the acceptance probabilities yields that the target distribution A^j^+n is 
left invariant by the Metropolis-Hastings dynamics. □ 



5 Application to linear network extraction 

The goal of this section is to apply the model presented in Section [2] to the extraction 
of a network of tracks in between crop fields from image data. The left-most panel in 



16 



Figure [6l obtained from the collection of publicly released SAR (Synthetic Aperture 
Radar) images at the NASA/JPL web site http://southport.jpl.nasa.gov shows 
an agricultural region in Ukraine. A pattern of fields separated by tracks is visible, 
broken by some hamlets. The image was previously analysed by Stoica et al. [25] by 
means of a Markov line segment process. 




Figure 6: Polygonal configuration (in yellow) overlaid on a SAR image of fields in rural 
Ukraine (top left panel), the corresponding edge map (top right panel) and a regular linear 
tessellation extracted from the Hough accumulation array (bottom row). 

Note that the tracks that run between adjacent fields show up in the image as 
whitish lines against the darker fields. Thus, a track is associated with a high image 
gradient. To find a suitable family T of straight lines, we therefore begin by computing 
the gradient of the data image after convolution with a radially symmetric Gaussian 
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kernel with standard deviation a = 3 to suppress noise. The right-most panel in Fig- 
ure [6] shows the gradient length thus obtained. We then compute the Hough transform 
[HI Ho] in an 80 x 80 accumulation array and select the 8 lines corresponding to the 
bins with the largest number of accumulated votes. This collection is augmented by 
lines parametrised by the largest local extrema in the accumulation array to yield the 
42 lines shown in the bottom panel of Figure |6l The resulting set is finite and contains 
no three lines that intersect in a common point, in other words, it is a regular linear 
tessellation. Regarding the choice of ay, a comparison of the data with the simulations 
shown in Figure [2] suggests ay = 1/2. Since the lines are selected on the basis of a high 
gradient value, there seems to be no particular reason to favour one line over another, 
so we may set the line activity to a constant value. 

To quantify how well a polygonal configuration fits the data, recall that an edge 
should be present when there is a large gradient, and absent when the gradient is small. 
This desirable property is captured by the Hamiltonian 

-HDil) = -(3 Yl [/(^) - ^(^)] ' (5) 

eG-B(7) 

where /(e) is the integrated absolute gradient flux along edge e, c(e) a threshold to 
discourage spurious edges, and /3 > 0. We take c(e) proportional to the number of 
segments along the edge with proportionality constant c > 0. For this choice, ([5]), 
being a sum of segment contributions, is local in nature, which is convenient from a 
computational perspective. 

To find an optimal polygonal configuration, we use simulated annealing applied to 
the Metropolis-Hastings algorithm of Section [3] for A^^^+-Hd with given by and 
Hd defined by and recolour rate 100. Starting at temperature 1//3 = 10, the inverse 
temperature parameter /3 is slowly increased to 100 according to a geometric cooling 
schedule. The result for c = 2, line activity 1/2 and k = 4 is shown in Figure [61 There 
are no false negatives; a few false positives occur near the hamlets and are connected to 
the track network. The precision of the line placement is clearly linked to the precision 
of the underlying regular linear tessellation. 

6 Summary and discussion 

In this paper, a new class of consistent random fields was introduced whose realisations 
are coloured mosaics with not necessarily convex polygonal tiles. The vertices of the 
polygonal tiles may have degrees two, three or four. The construction, inspired by the 
Arak-Clifford-Surgailis polygonal Markov fields in the continuum [21 [3], [1] , extends our 
previous construction ^23] of consistent polygonal random fields with tile vertices of de- 
gree two only. The latter case is substantially simpler due to the fact that interactions 
are restricted to a hard core constraint only, and, moreover, for a given realisation, 
there are only two equally likely admissible colourings. We developed a dynamic rep- 
resentation for consistent multi-colour polygonal fields, which was used to prove the 
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basic properties of the model including consistency and an explicit expression for the 
partition function. Local and spatial Markov properties were also considered. The dy- 
namic representation provided the foundation on which to build Metropolis-Hastings 
style samplers for Gibbsian modifications of these fields. Finally, we applied the model 
to the detection of linear networks in rural scenes. A modification of our models would 
consist in ascribing activity parameters to segments rather than lines, cf. |T7]. A dis- 
advantage seems to be that the dynamic representation would depend on the direction 
of time, in other words, the model would be anisotropic. 

To conclude, we should stress that, although the model was inspired by those of 
Arak et al. there are striking differences inherent to the discrete set-up. Notably, in 
the continuum collinear edges are not allowed, in the discrete set-up they are. Indeed, 
if one were to forbid collinear edges, this would lead to a forbidden line whose influence 
would be felt at arbitrarily large distance from its single edge, hence ruling out any 
meaningful Markovianity. This is not true in the continuum as the dynamic represen- 
tation there ensures that collinear edges occur with probability zero. As a consequence, 
our consistent random fields are not Arak-Surgailis fields conditional on having their 
edges along the lines of T. More fruitfully, our models provide a graph-theoretical in- 
terpretation of mutually compatible Gibbs random fields that inspires novel simulation 
algorithms as an alternative to the usual local tile updating schemes [29]. 
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